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ABSTRACT 

The design of high order control systems by parameter 
optimisation usually involves a great deal of computation. 
Round-off errors in machine computation may lead to inac- 

. v‘ 

curate results that in many cases may prevent the continu- 
ation of the optimization procedure.;""' Some of the numerical 
operations involved in a currently used optimization tech- 
nique are discussed and analyzed with special attention to 
the numerical accuracy. Alternative methods are suggested 
for more accurate and cost effective solutions. 
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CHAPTER 1 


BACKGROUND AND PROBLEM IDENTIFICATION 


1 . 1 The Design of a Control System by Parameter Optimization 

The design problem that raised the need for this work 
is the one of finding the optimal set of prespecified para- 
meters of a fixed configuration control. system. A criterion 
for the design can be derived by some evaluation of the sys- 
tem response to some standard input. If the system's trans- 
ient response to a step input is weighted by some criterion 
and integrated along all time, a performance ind^x can be 
defined, the minimization of which gives the optimal design. 

In the case of a linear system, whose state vector sat- 
isfies the equation: 


x ( t) = Ax(t) x (0) = Xq (1.1) 

"'pii-’-,. 

such a performance index can be given the quadratic form: 


T 

PI = x A w x n + 

— u — u 


fi 


T 

Q x dt 


( 1 . 2 ) 


where x is the transient state vector, is the initial 
transient state vector and W and Q are positive semi- defin- 
ite weighting matrices. 

The integral expression in (1.2) can be presented more 
conveniently by: 
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(1.3) 


/ OO 


x Qx dt = tr [Q* X] 


where tr denotes the trace of the matrix, and where: 


X = 


T 

x x dt 


(1.4) 


Using equation (1.1): 

r3 T> « rp * m TT T 

dt ^ “xx 1 + X X = X X A + AX X (1.5) 

and integrating both sides along all time for a time invari- 
ant system gives: 

; v 1 - ” 

x x T (t=°°) - x x T (t=0) = XA T + AX 

For a stable linear system, the transient response vanishes 
as t and we get that the state matrix X must satisfy 

the equation: 


AX + XA ~ -X 


0 


( 1 . 6 ) 


where : 


X 0 = -0 2Eq‘ 


(1.7) 


1 . 2 The Model Performance Index 

The Model Performance Index was first suggested by 
Rediess “[ll] and was improved by Palsson [10] for applica- 
tion in an automatic optimization algorithm. The theory is 
developed in these two references and will not be described 
here. A brief outline of some results, applicable in an 
optimization procedure is necessary, however, for the under- 
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standing of some computation operations analysed in the 
following chapter. 

It is desired to approximate the step input response of 
a system whose transfer function is : 


, m , 
b s + 

G (s) » — 

s n 


+b^s + b Q 


2 ^ 

s" + a^_^s . +a^s+ 


( 1 . 8 ) 


to that of a model whose transfer function is: 


3 v s + +3 s + S n 

G m (s) = r - r i 2 

s + _ j_ s +. . . . +a^s+ag 


( 1 . 9 ) 


The model is required to satisfy the condition: 
l - k <_ n - m 

A new transfer function that has no zeros, and whose denomin- 
ator is the cascaded system denominator and model numerator, 
is formed: 


G(s) = 


a 0^0 


(e^s k + +3-,s+3 n ) (s n +a n ^ 1 s n " 1 -f. . . +a n s+a n ) 


with the initial conditions: 

b. 


x 


_0 

x o a o 


x , . . , x - 0 
(x+X) 0 


1 <i < n-m 


i-1 


x (i+l) = b n-i - Z a n-i ' J x ' ' 


i+j (j+1) 


j=n-m 


1 w 0 

( 1 . 10 ) 


n-mcicn 


( 1 . 11 ) 
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The Model Performance Index is defined as: 


where : 



x^Qx dt 


Q « 


a. a 

• a o“o 


a. 12) 


(1.13) 


d is the model characteristic coefficients vector presented 
in the n'th dimension state space : 

d — [otQ , ■ . . a ^ 1 r -0 1 0 * • • 0 ] (1.14) 

x in the equation (1.12) is the transient state Vector of the 
system represented by equation (1.10) in response to the in- 
itial conditions (1 . 11) . By use of expressions derived in the 
previous section, the performance index can be expressed as: 


PI = tr [Q *X] 


X must satisfy the matrix equation 


AX + XA - -X 


0 


where : 


X 0 = 


(1.15) 


( 1 - 6 ) 


(1-7) 


1 . 3 The Solution of the State Matrix X 

The solution of equation (1.6) for the matrix X is a 
major concern of this thesis. Palsson [qo] has shown that 
when the system matrix A is written in the phase variable 
form: 
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0 10 0 0 

0 0 1 0 0 



(1.16) 


the first column of the state matrix X, x^,can be computed 
by: 


n 


x = -E ^ / E. , x, 
— 1 n l-l — ( 


i=l 


n-i 


(1.17) 


where x^ is the (n-i) 1 th column of the matrix X~ and the 
n-i u 

matrices E^ are defined by the relationship: 


with 


E . = -AE . + a , I 

i l-l n-i 


E o = 1 


l<i<n 


(1.18) 


The remaining columns of the matrix X are computed by the 
relationship : 


x. = -Ax. + x rt 
—l -l-l — 0 . 

i 


(1.19) 


The computation of these expressions is discussed in detail 
in Chapter 2. 


1 . 4 The System Transfer Function Coefficients 

A method for obtaining the transfer function of a com- 
licated (multi-loop with a multi-input multi-output control- 
led member) system given by its components was developed by 
Whitaker [13], Griffen [4] . and Beyers [1]. It consists of 
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the following generalized steps: 

1. Obtaining the system state equations in the form: 

x = Fx + Gu 

y “ Hx + Ju (1.20) 

where x is the state vector, u is the input vector, 
is the output vector and F, G, H, and J are properly 
dimensional matrices. 

2. Obtaining the transfer function matrix by: 

%_ = [H[sI~F] _1 G + J] u (1.21) 

This method is used to obtain the system coefficients which 
make up the matrices A and X Q in equation (1.6). Although 
the accuracy of the computation at this stage is not analysed 
in detail in this work, it should be pointed out that the 
large number of arithmetic operations done here especially 
for high order systems, has been found to be a source of 
significant numerical errors. The accuracy of the computed 
performance index cannot be any better than that introduced 
in the matrices A and 

The error in the system coefficients may affect the min- 
imum search procedure mainly in two ways : 

1. When the performance index first increases in a step, 
it is normally assumed that the minimum point lies 
somewhere along this step. In the presence of signifi- 
cant inaccuracies such an interpretation of the perform- 
ance index increase may be false. Taking a shorter step 
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in this case may even give further increase of the per- 
formance index (this problem is later referred to as 
the "performance index increase in half a step"). 

2. In the current optimization program used at the 
Measurement Systems Laboratory at M.I.T. and applying 
the methods discussed in this chapter, the gradient 
method is used for the search of the optimum. The 
derivatives of the system coefficients with respect to 
the design parameters are used to produce the gradient 
at each step. These derivatives, which are computed by 
incrementation of the parameters and division of the 
corresponding changes in the coefficients by the incre- 
mented parameters values, are very sensitive to the 
accuracy in the computed parameters. In other words, 
if the error in the computed coefficients is of the 
order of the expected incremental change in the coeffi- 
cients , then the resulting derivatives and gradient 
values would be completely incorrect. 

These two problems are considered again in Chapter 4 where 
an alternative method is discussed. 

1 . 5 Scaling the System Coefficients 

A large number of operations involving matrices and 
vectors take place in the solution of equation (1.6). In some 
cases significant inaccuracies may occur when elements appear 
in certain arrangements of magnitude in matrices and vectors.. 
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These cases are discussed in Chapter 2. Scaling these 
matrices and vectors so as to rearrange the relative magni- 
tudes of the elements may improve the results significantly. 
Since all the elements in equation (1.6) consist of system 
coefficients, the desired scaling of matrices and vectors 
can be obtained in most cases by scaling the system coeffi- 
cients. 

One way of scaling the coefficients is changing the time 
scale. Since the units of the coefficients are powers of 
time, changing the time scale can change the magnitude of 
the coefficients in the desirable way. This is illustrated 
in the following: 

Characteristic coefficients before scaling: 

a 0 / .... a n _ 1 ,1 


Characteristic coefficients after time scaling: 


-n ^n-1 ^ -i 

f a Q , f ....f a n _]_/l 


where f is the scaling factor. By selecting a value for f, 
the coefficients can be ordered in an ascending or a descend- 
ing order of magnitudes as desired. This type of scaling 
will be referred to later as time scaling. 

Another way of scaling the system coefficients is add- 
ing a pair of a pole and a zero of the same value (or at the 
same location in the complex plane) which, of course, does 
not change the characteristics of the control system itself. 
The relations between the characteristic coefficients 
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and the poles are as follows : 

a Q = product of all poles 
a^ = sum of all possible products of n-1 
out of n poles 

a 2 = sum of all possible products of n-2 
* out of n poles 

« 

a^ ^=sum of all poles 

It follows that adding a pole-zero pair near the origin, for 
example, would decrease a Q , not affecting significantly, 

since : 

a 0 ■ Pl-?2 P n ' £ 

and 

a -i " p_+p,_ + p +e 

n-1 

where P]_««*P n are the original system poles, and £ is the 
additional scaling pole (assuming that e is not of the mag- 
nitude order of the other poles) . There is freedom in the 
location of the pole-zero pair and the corresponding effect 
on the coefficients. This type of scaling will be referred to 
later as pole-zero scaling. There are some operations in the 
solution of equation (1.6) (as discussed in sections 2.2 
and 2.3) where time scaling of the coefficients cannot im- 
prove the accuracy, but pole-zero scaling may be used effec- 
tively for this purpose. The disadvantage of the pole-zero 
scaling is that it increases the order of the system. The 
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trade-off between the positive effect of scaling and the neg- 
ative effect of increased order on the accuracy must then be 
considered. 

1 . 6 Severe Inaccuracies Occurance 

There were some typical cases in which severe inaccura- 
cies occurred, preventing the completion of the optimum search 
procedure. It is impossible, however, to point out one gen- 
eral reason for the occurance of the problem. Severe inac- 
curacies result in phenomena like negative performance index 
values, the ’’performance index increase in half a step", 
large values of performance index when small ones are expect- 
ed, and other results that just do not make sense. 

The problem has occurred in cases where the complex fre- 
quency range of the system modes of response was large, or 
in other words, the system poles were far apart. One such 
example is the presence of the phugoid mode in the equations 
of a controlled aircraft* This mode has relatively low frequen 
cy and small magnitude. When the "short period approximation 
was exercised and the phugoid mode was eliminated, the severe 
inaccuracies in the computation disappeared. The addition 
of a remote pole to a system may also cause the appearance 
of severe inaccurracies in the computation. In such cases 
the problem can be avoided by elimination of low residue 
mqdes. But this is not always the case. In many situations 

low residue modes are not easy to identify and to isolate, 
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especially in the process of optimization, and it is not 
always clear that such elimination would solve the inaccura- 
cy problem; Thus, rather than treating the different cases 
separately, the approach to the problem was by investigation 
of the numerical operations trying to identify the potential 
sources of inaccuracies and to suggest adequate solutions. 
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CHAPTER 2 


THE NUMERICAL COMPUTATION 


2 . 1 Introduction 

The method proposed by Palsson for the solution of the 
matrix equation (1.6) was presented in section 1.3. In this 
chapter each step of the solution procedure is discussed in 
terms of the explicit algebraic expressions and the arithme- 
tic operations, with special attention to the numerical ac- 
curacy obtained in these operations. Some parts of the com- 
putation are discussed only for the purpose of shedding some 
light on the various stages of the procedure, so that . 
sources of significant inaccuracies that may occur can be 
spotted more easily. 

2.2 The Computation of the Matrices l<i<n 

The matrices are obtained by the relationship: 

E^ = -AE i _ 1 + A n _^I l<i<n (1.18) 


with 


= I 


This part of the computation can be best investigated by con- 
sidering the algebraic expressions of the elements of these 
matrices. The first few matrices are given in Appendix A so 
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that the nature of these and of higher order E^ matrices can 
be understood. 

From the explicit expression of the E^ matrices, it can 
be learned that the elements at lower positions in the col- 
umns contain higher powers of the system characteristic 
coefficients. The arrangement of the coefficients in the 
matrices elements may affect the accuracy of the computation 
in two ways : 

1. inaccurate computation of the matrices elements; 

2. ill conditioning of the matrix E n with respect to 
its inversion in a later stage of the procedure. 


To understand the first possibility it should be realiz- 
ed that in the computation of low position elements of the 
E^ matrices , small and large numbers are normally added to- 
gether with considerable loss of low position digits. 

Notice y for example, that the element: 


E 5 (5,5) = 2a n _ 5 + 2a n _ 3 a n _ 2 - 2 ^^ + 2a^_ 2 ^_ 


2 3 5 

+ 6 a 0 a , • — 6 a ~ a 1 + 2a , 
n-3 n-1 n-2 n-1 n-1 


would be normally dominated by the higher powers of a n ^ , un- 
less some scaling of the coefficients is done. Notice that 
time scaling of the system coefficients would have no effect 
on the accuracy of this computation, since such scaling is 
equivilant to multiplying the matrix element by a factor, 
leaving the relative magnitudes of the addents unchanged. 
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Using a scaling factor f gives, for example: 
E 


5 (5,5) = 2 ( f 5 a 5 ) + 2(f 3 a n _ 3 f 2 q n _ 2 ) -2 
scaled 


+ 2 <<f 2 a n _ 2 ) 2 fa n _ 1 ) + 6(f 3 a n _ 3 (fa n _ 1 > 2 ) 
- 6(f 2 a n . 2 (fa n _ 1 ) 3 ) + 2(fa n _ 1 ) 5 


= f J E 5 (5,5) 

unsealed 


Pole-zero scaling can be used effectively to control the 
accuracy of this computation. Such scaling, as was dis- 
cussed in "Section 1.5, can be used to magnify the value of 

a A and other low order coefficients with respect to a n 
0 n-1 

without changing the value of a n _^ significantly. Scaling 

would be done to make the addents close enough in magnitude 

so that loss of low position digits is minimized. 

The inversion of the matrix E is discussed later in 

n 

this chapter. It is pointed out that ill conditioning of a 

,Y' 

matrix - with respect to its inversion occurs when the ele- 
ments on the diagonal are significantly smaller than ele- 
ments to their right in rows. It is evident from the ma- 
trices (consider E^,for example) that such situations can be 
controlled by either time scaling or pole-zero scaling of 
the system coefficients. A situation where elements on the 
right of the rows are larger than those on the left can be 
reversed by such scaling. 

One may argue that changing the order of the operations 
in the computation of the matrices so as to sum small nuim 
bers first and then add the larger numbers, minimizing loss 
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of digits, would be useful. It can be verified, however, 
that in general these operations are indeed done in such 
order, as each E^ matrix is computed from the previous one. 


n 


2.3 The Computation of X) E j_-i^o 

i"l n-i 


n 


- 0 . 


in the expression X) E i-l r — 0 1 th 


n-i i-1 "n-i 

column of the initial conditions matrix: 


x 0 = 2<£o 


X 0lX 

’Vo, 


x o x o 

U 1 u 2 


X 0. x o 

1 n 


x o x o. x o x o 9 

n 1 n 2 


x o x o 

n n 


( 2 . 1 ) 


The initial conditions are computed by equations (1.11) as 
explicit functions of the system coefficients. Consider the 
case: n=8 , m~5, t-3 , then: 


x 


X 


X 


X 


X 


0 

0 

0 

0 

0 


1 

2 

3 

4 

5 



0 

0 

b 

b 


5 

4 
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X 

- b-, - a^b r - a_b „ 
3 6 5 7 4 


X 

= b 2 - a 5 b 5 ‘ a 6 b 4 ' a 7 b 3 + a 7 a 6 b 5 

+ a«a„b. 
7 7 4 

CO 

o 

X 

b l a 4 b 5 " a 5 b 4 “ a 6 b 3 + a 6 a 6 b 5 

+ a - a^.b . 
6 7 4 


a 7 b 2 + a 7 a 5 b 5 + a 7 a 6 b 4 + a 7 a 7 b 3 

— a„a_a^b 
7 7 6 


- a^a^a^b , 
7 7 7 4 


( 2 . 2 ) 


The higher order initial conditions involve higher order 
products of the system's coefficients, and normally they are 
significantly larger than the lower order ones. Computation 
of high order initial conditions involves a great deal of 
multiplications and additions , and single precision is in- 
adequate if accurate results are desired. Since each one of 
the initial conditions depends linearly on one b coefficient, 
it can be argued that inaccurate computation would result in 
initial conditions corresponding to another set of b coeffi- 
cients or , in other words, to another system. This may 
cause problems like the "performance index increase in half 
a step" , and inaccurate computation of the gradient of the 
performance index with respect to the design parameters, 
that were discussed in Chapter 1. The accuracy of x Q is of 
special importance. Since in the optimization algorithm a 

system static sensitivity of 1 is assumed, the value of x„ 

°1 

must be 1. Even a slight error in this value would result 
in a steady state difference between the model's and the 
system's responses, which, at least theoretically, should 
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produce an infinite performance index. This point was not 


considered in the original program where the value of x 


0 


was computed in 3 steps, which in single precision floating 
point arithmetic resulted in 2 incorrect digits out of 8. 

This error has lead to "awkward" values of performance index. 
When the value of x n , was fixed at -1 (instead of the unnec- 
essary computation) sufficiently accurate results were ob- 
tained. 

The error in x n had a major effect on the computed per- 

1 

formance index in the case of one system, but had no signifi- 
cant effect in the case of another system of the same order. 

This is explained by the "weighting" of the elements x^ x^ 

u i u j 

of the Xq matrix by the E. matrices, first in the sum 

n 1 -1 n 

Ee. , x n and then in the product E„ E E. x. 
i 1 n i i — 1 n — i 

As was shown in 2.2 the lower order E^ matrices contain many 

zero elements. It is easy to realize that in the vector 

.. r': ' 

n 

E E i-1 —0 first element is mostly affected by pro- 

i=l ~ n-i 

ducts of x n . Elements at the lower positions in the vector 
1 . 

are less affected by x^_ . This effect depends, of course, on 
the relative magnitudes of the elements of e. , or the weight- 

l 

ing imposed by these matrices. Then, in the product 
_ x n 

E E E -i i x n T the contribution of the elements of the 
n • i—i — u j 


i=l 


n-i 


n 


-1 


vector E E. .. depends on the weighting by E ^.When 

i-l 1-1 n-i 

— 1 

the first elements in the first row of E are considerably 


25 



larger than other elements, the error in is weighted most 
heavily in the first element of the product vector x^. This 
indeed was the case when the problem was first encountered. 
Fixing the value of x^ at -1 resulted in a significant 
change in X(l,l) (x-^(l)), the other elements remaining prac- 
tically unchanged, which in turn gave the desired correction 
in the computed performance index. 

The algebraic expressions of the elements of -^Xq 

n-i 

are extensive and giving them here in detail seems unneces- 
sary (besides, they can be easily obtained from the given 
and Xq matrices) . It should be noted, however, that this 
part of the computation involves a great deal of arithmetic 
operations and may be a source of severe inaccuracies. If 
these inaccuracies are found to influence the computation 
significantly, rearranging the operations so that small and 
large numbers are summed in an optimal order (to minimize 

, i'V . 

loss of low-position digits) should improve the results. If 

this is done in the computation of the product vectors 

E i-1— 0 only, the additional storage required is negligable 

1 n-i 

(n words) , but the additional check operations may be time 
consuming. The alternative of extended precision computation 
has the disadvantage of having to use a different program- 
ming language (see 2.4) if precision higher than double is 
desired. The cost of interfacing languages in a program is 
considerably high, and the trade-off between the different 

possibilities must be considered. 
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Considering the corresponding expressions in the mat- 
rices E^ and Xq, it can be verified that time scaling of the 
system's coefficients is eqUivilant to multiplying the re- 
sults by a factor, and thus, has no effect on the accuracy 
in the computation of the vectors or of the sum 

2 E i-1 ~o (see a verification of a similar argument in 
i=l n-i 

2.2). Pole-zero scaling can be used effectively to arrange 
this computation so that small numbers are summed first and 
then added to the larger numbers. This can be done by 
placing the pole-zero pair on the left half of the real 
axis far from the origin which would make appreciably 
smaller than a Notice, however', that such scaling may 

have an undesirable effect on the E matrix, increasing the 
elements to the right of the diagonal in the rows. As is 
discussed in section 2.4 this may "ill condition" the matrix 
with respect to inversion. 


2 . 4 The Computation of and Matrix Inversion 

is the first column of the state matrix X- 
isfies the equation: 


n 

E x, = E . , x n 

n—1 l-l —0 

i=l n-i 


It sat- 


(2.3) 


or 


^ ■ '- 1 & ^ -°»-i 


(1.17) 


This section is mostly concerned with the numerical errors 
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arising in the inversion of the matrix for the solution 
of X . 

A popular method for machine inversion of matrices is 
the Gauss-Jordan method. This method is described in de- 
tail in references [5] , [6] , [7] , and [12] and will not be 

presented here. It is , however, the matrix inversion 
method referred to in the following discussion. 

A significant loss of accuracy may occur when elements 
on the diagonal of the matrix to be inverted are consider- 
ably smaller than elements to their right in the rows. 

Such situations may be avoided by reordering the rows of the 
matrix so as to place the larger elements on the diagonal. 

Such reordering may require extensive check and bookkeeping 
operations for high order matrix inversion. In the case of 
the matrix such a solution is inadequate in general, 
since elements in all the rows appear in a similar arrange- 
ment -bf relative magnitudes (see Appendix A). 

The usefulness of scaling a matrix by a linear transform- 
ation before scaling is questionable. The scaling and rescal- 
ing transformations may cause additional errors. Multiplying 
the matrix elements by powers of 10 when floating point arith- 
metic is used may result in a matrix that is isomorphic to 
the original. Indeed, attempts to use such scaling by a trans- 
formation matrix did not improve the inversion accuracy. 

Scaling the system characteristic coefficients is, as 
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discussed in section 2.2, a useful tool for initially arrang- 
ing the elements of the E matrix in a desirable order of 

n 

magnitudes. This is illustrated in an example case in Chap- 
ter 3 . 

The situation of diagonal elements that are considerably 
smaller than elements to their right may arise at a later 
stage of the inversion procedure. How the pre-conditioning 
of the matrix should be done to avoid such cases is not 
easy to determine, since the relation between the elements 
of the resulting matrix and the elements of the original ma- 
trix is not a simple one. One way to handle such situations 
is suggested by Hellerman [5, p.57J.' It consists of using 
the largest element in the corresponding column rather than 
the diagonal as a pivot element, while recording the place 
of this element* Then, since the rows and the columns of 
the inverse matrix appear scrambled, rearranging them in 
the right order (the Gauss -Jordan method is referred to as 
the inversion method) . 

A number of measures, and criteria for the inaccuracy in 
matrix inversion and linear equations solution procedures 
can be found in the literature. They all fail to give an 
exact measure of the inaccuracy (to indicate how many digits 
in a result number are incorrect) and, at best, give an indi- 
cation that some significant inaccuracy has occured in the 
computation. One such check often used is obtained from the 
elements of the difference matrix: 

" A* A 
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In our case of interest this check is hardly satisfactory. 

An investigation has shown that when the inverse of con- 
tains elements considerably different in order of magnitude, 

the small elements are computed with large errors. In the 

-1 n 

computation of the product E ^ E._-. x n the contribu- 

1=1 n-i 

tion of these elements is negligible, and even relatively 
large errors in them would have no effect on the final re- 
sult. Yet, when the product E^ is computed for the 

check, these elements are multiplied by large elements of 
the matrix E , and the errors in them may produce appreci- 
ably large elements in the check matrix D. On the other 
hand, errors in the large elements of E \ which can affect 
the computation of significantly, may not be detected by 
the check at all. Thus, such checks can be used to obtain 
a general idea about the accuracy of the computation but 
not as f a measure of it, especially not as a basis for com- 

t . t'V. '.- r 

parison between numerical results. 

Absolute large values of the matrix elements (not only 
their relative magnitude with respect to each other) have 
been found to affect the accuracy of the matrix inversion. 
Severe inaccuracies have occured in some practical cases 
when the matrix determinant was very large. When the ele- 
ments of the matrix were divided by a certain factor before 
and after the inversion, accurate inversion was obtained in 
these cases. 

A different approach to the solution of the vector x^ 
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in the matrix equation: 


^ ■ £ 


E . n X. 
i-l — 0. 


n-i 


(2.3) 


is usually referred to as the indirect or the iterative ap- 
proach. The Gauss -Seidel method is a popular iterative meth- 
od used for machine solution of systems of linear equations. 

It is described along with some other such methods in refer- 
ences [6] and [7], where it is also pointed out that these 
methods may or may not converge on the solution depending 
on the numerical values involved, even in the absence of 
round-off errors . 

If a sufficiently accurate computation of x^ cannot be 
obtained by scaling the matrix, an extended precesion ver- 
sion of the matrix inversion algorithm may be used. Since 
IBM/370 FORTRAN does not have a built” in extended precision 
facility, the MINV program from the IBM PL/1 Scientific 
Subroutine Package [8] has been modified to invert a matrix 

. V.’ 

in quadruple precision. This program can be called by the 
FORTRAN optimization program whenever an extended precision 
inversion is desired. The modified program and the inter- 
face control cards are given in Appendix B. 


2 . 5 The Remaining Columns of the Matrix X 

After the first column of the state matrix X is comput- 
ed, the remaining columns are computed by equation (1.19) : 


*i - -a*!.]. + 5o. 

1 
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e . g. 

0 

-1 

0 

0 ... 

r 

o 


vo; 

-2 

0 

4 

1 

• 

0 

-1 

0 . . . 

0 

Si + 

« 

• 


a o 

a i 

• • • 
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x o x o , 

n 1 
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Note that except for the last element in each column all 
the elements are obtained by shifting the elements of the 
previous column upwards with an inverted sign and adding an 
element from the previously computed Xq matrix. Both oper- 
ations do not produce any numerical error. Figure 2.1 the 
elements of an 8 ' th order state matrix X that are computed 
by shifting the elements of the first column and adding an 
element from the Xq matrix are denoted by x, and elements 
that involve multiplication of the system's coefficients by 
a previous column are denoted by *. The computation of the 
* elements in this stage from previously computed values 
may add. to their inaccuracy. 

xxxxjxxxx* 

xxxxixxx** 

XXX x] XX* * * 

X _x__x_ X! X * * * * 
xxxx***** 

XXX*** * * * 

XX******* 

x******** 

_J 

Figure 2.1 - Additional Inaccuracy in the Computation 
of the Remaining Columns of an 8'th order X Mateix 


In the case of the Model Performance Index only the elements 

32 



that are both in the first £+1 rows and in the first £+1 
columns ( £ being the model's order) take part in the com- 
putation of the performance index (e.g. in the case of an 
8 ’ th order system and a 3 ' rd order model only the elements in 
the 4x4 square in Figure 2.1 affect the result). It fol- 
lows that for most practical cases this part of the compu- 
tation does not introduce any numerical inaccuracy at all. 

2.6 The Computation of the Performance Index 

The performance index is computed by equation (1.15) : 

PI = tr [Q-X] 

where Q is the weighting matrix in the performance index and 
X is the state matrix whose computation was discussed before. 
Except for the case of very high order systems and when the 
correct value of PI is close to zero, this part of the com- 
putation should not introduce additional severe inaccuracy. 

If the expected value of PI is small, then accurate low-po- 
sition digits must be obtained, and single precision may not 
be sufficient. Also when the order of the system is high, 
the large number of multiplications and additions may result 
in significant loss of digits. Note that in the case of the 
Model Performance Index the elements in the matrix Q that 
are both not in the first l+l rows and not in the first £+1 
columns are all zeros, and thus, the amount of computation 
in this stage depends only on the order of the model. 
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CHAPTER 3 


AN EXAMPLE CASS 


3-1 Introduction 

Because of the nature of the problem the numerical ex- 
ample cannot be of a simple, low order system. A real in- 
accuracy problem encountered in a control system design is 
discussed, in order to demonstrate some of the inaccuracy 
problems and their solutions as described in the previous 
chapters . 

3 • 2 The Example System and the Inaccuracy Problem 

The example control system is given by its mathematical 
block diagram in Figure 3.1. The numerical values are given 
in Figure 3.2(a) in the format of the System Description Pro 
gram [12] . The rows of the SYS and the SIG matrices corres- 
pond to the blocks in the block diagram. The last two rows 
of these matrices represent the model zeros that are cascad- 
ed to the system to use the Model Performance Index. The 
controlled member state matrices and the design parameters 
location and initial values are also given. The computed 
values of the over-all transfer function characteristic co- 
efficients (ACOF(I)), the numerator coefficients (BCOF(I)) 
and the initial condition vector are listed in Figure 3.2(b) 
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The transfer function coefficients are time scaled by a fact- 
or of two. 

Significantly inaccurate computation was experienced in 

attempting to optimize this system. The numerical values of 

the matrix E , the determinant of this matrix and the state 
n 

matrix X corresponding to the initial values of the design 
parameters are given in Figure 3.3 for comparison with other 
cases. The computed performance index in this case was 0.052 
while the correct value is 0.26, as was verified by integra- 
tion methods. Notice that the computed value of the first 
element of the matrix X is 0.55. The correct value of this 
element was found to be 0.75. The other elements of the 
matrix X are sufficiently accurate. 

3 . 3 The Practical Solution of the Inaccuracy Problem 

In order to obtain more accurate results in the computa- 
tion of the example case the value of the first initial con- 
dition was fixed at -1, and since the value of the E matrix 

n 

determinant was very large, all the matr,ix elements were de- 
5 

vided by 10 before the matrix inversion. The resulting 
and X matrices are given in Figure 3.4. The given elements 
of the matrix are before reduction for comparison with 
other cases and the determinant was computed after reduction. 
As can be seen in the figure, the computed value of the first 
element in the matrix X is now 0.75 while the other elements 
are considerably closer to the ones before the correction. 
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The resulting performance index value was 0.26, differing by 
less than one percent from the value obtained by integration. 
Similar agreement between the performance index values com- 
puted by integration and the ones computed by the method 
discussed in this work was experienced in the following steps 
of the optimization . It is interesting to note that the 
elements at lower positions in the columns (or the rows) of 
the X matrix are closer to the correct values , which is in 
agreement with the error weighting effect, discussed in sec- 
tion 2.3 

The correction of the initial condition has, of course, 

no effect on the matrix E . This correction affects the sum 
n n 

/] E. -.x~ while the reduction of the matrix E elements 
i-l—O . n 

i—l n - 1 

affects the matrix inversion, so that the first column of the 

-1 n 

matrix X which is equal to E Y] E. -,x n is affected by 

n -i , i— l— u 

i=l ' n- i 

both corrections. Indeed, it was found that both corrections 

were necessary in this case. When only x n was corrected 

1 

and the reduction of the E matrix elements was not exercis- 

n 

ed, negative performance index values were obtained at later 
steps of the optimization. 

3 . 4 Scaling the Transfer Function Coefficients 

As was discussed in the previous chapters scaling the 
transfer function coefficients may be useful in conditioning 
matrices and vectors for more accurate computation. Consid- 
er the case where a pole -300 was added to the example system 
presented in section 3.2. The optimization program current- 
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ly used at the Measurement Systems Laboratory at M.I.T. re- 
duces the system coefficeints by time scaling to prevent large 

values of the E matrix elements which are undesirable as 

n 

discussed previously. The scaling that was done in the re- 
mote pole case is shown in Figure 3.5. This scaling has re- 
sulted in an unfavorable conditioning of the matrix E^ where 
some of the diagonal elements are smaller than elements to 
their right, as shown in Figure 3.6. 

When the scaling factor was inverted as shown in Figure 

3.7 and the increase of the E matrix elements was handled 

n 

..-■•k 15 , . 

by dividing them by 10 , the resulting condition of the 

matrix was a desirable one, as the elements on the diagonal 
dominate in magnitude the elements on their right. This is 
shown in Figure 3.8. The resulting error in the computed 
performance index was in this case 30% of its value. 

It is not argued that changing the condition of the ma- 
trix E would usually improve the accuracy of the computation. 
In sbitie cases it was found to have no effect at all. The 
possibility to control the condition of the matrix by such 
scaling should be noted and used when necessary. However the 
problem of large values of the E n matrix elements should be 
treated directly by dividing these elements by a factor and 
not by scaling. Scaling should be used on a selective basis 
to condition matrices and vectors for specific operations 
like the matrix inversion. The establishment of more specific 
criteria for the usefulness of such scaling in different cases 
would provide the control system designer means to obtain more 

accurate computation in many practical cases. 
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Figure 3.1-A Mathematical Block Diagram and Data Input Definition of the 

Example Control System 
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Figure 3 . 2 (a) -Input Data for the Example System in the System Description Program Format 
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Figure 3 . 2 (b ) -Computed and Scaled Transfer Function Coefficients and Initial 

Condition Vector 
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Figure 3.3-The Matrices E^ and X Corresponding to Initial Parameter Values of 

the Example System 
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Figure 3 . 4~Corrected Values of the Matrices E n and X of the Example Case 
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Figure 3.5'Scaling the System Coefficients of the Remote Pole Case 
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CHAPTER 4 


ALTERNATIVE METHODS FOR APPLICATION IN A PARAMETER 
OPTIMIZATION PROGRAM 


4.1 Introduction 

The search algorithm of a parameter optimization program 
must utilize other routines that execute some specific parts 
of the computation, like the system mathematical representa- 
tion and the solution of equations like the matrix equation 
(1.6). The methods presented in this chapter are an attempt 
to take more accurate and cost effective approaches to the 
problems of the matrix equation solution and the system re- 
presentation. 

4 . 2 A Minimal Method for the Solution of the Matrix Equation: 

AX + XA T = ~X Q 

The method for the solution of equation (1.6) that was 
suggested by Palsson [10] and was described and analyzed in 
sections 1.3 and 2.1 through 2.5 is extensive and for the 
most part inefficient. The number of arithmetic operations 
in this method is very large, which makes the solution of the 
equation for many practical design problemns both expensive 
and inaccurate. 

The method proposed in this section minimizes the num- 
ber of arithmetic operations needed to solve the equation. 
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It takes advantage of the fact that the matrices X and Xq 
are symmetric and the matrix A is in phase variable form. 

Since the matrix X is a symmetric n ' th order matrix {n 
being the order of the system) , there are — ■ - unknowns. 
Consider the case of a 4 ’ th order system. The matrix equation 
in this case is: 



X 11 X 12 X 13 x 14 
X 12 X 22 X 23 X 24 
x 13 X 23 X 33 X 34 
X 14 X 24 X 34 X 44 



where the system coefficients and the elements x Q of the 
matrix are known, and the elements of the matrix X 

are unknown . 


The following 10 independent equations for the solution 
of the 10 unknowns are readily obtained from the matrix 
equation : 
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2x 12 



= 

X 22 + 

X 13 



x 23 + 

X 14 


= 

X 24 

a 0 X ll ” a l X 12 " a 2 X 13 

a 3 X 14 

= 

2x 23 



— 

x 33 + 

X 24 


= 

X 34 

a 0 X 12 ” a l X 22 " a 2 X 23 

a 3 X 24 

— 

2x 34 



= 

X 44 

a 0 X 13 “ a l x 23 “ a 2 x 33 

a 4 X 34 

— 

-2a Q x 

14 “ 2 a l X 24 “ 2a 2 X 34 

2a 3 X 44 



“X 


-x 


-X 


-X 


-X 


-X 


-X 


-X 


-x 


-X 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


11 

12 

13 

14 
22 

23 

24 

33 

34 
44 


The equations are slightly modified and reordered in groups: 


group (a) 



1 

X 12 

2 X 0 

11 

1 

X 23 

*“ ““X a 

2 0 2 2 


1 

X 3 4 

" 2 X 0 
z u 33 


group (b) 
group (c) 


group (d) 


X 14 


-X 23 

x o 

u 13 


X 13 

= 

" X 22 

x o 

U 12 


X 2 4 

— 

~ x 33 " 

X 0 

u 23 


X 24 

- 

a o x n 

- a x x 

12 

X 2 3 

- 

a 0 X 12 

- a^x 

22 

X 44 

- 

a 0 X 13 

- a-^x 

23 

-2a 0 

X 14 

" 2a l 

X 24 

2 


2 X 13 

a 3 X 14 

— *0 
U 14 

2 X 23 

a 3 X 24 

=_x o 

u 24 

2 X 33 

a 3 X 34 



34 

2 X 3 4 " 2a 3 X 44 = ~ X 0, , 

4 4 


The original set of 10 independent equations is easily re- 
duced to a set of 4 independent equations in the following way: 
The unknowns in group (a) are now known. The unknowns in group 
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(b) ( here only x^) are directly obtained by substitution 
of solved unknowns of group (a) . The unknowns on the left 
hand side of group (c) and all the unknowns that have been 


solved are substituted into the equations of group (d) to 
give the reduced set of 4 equations for the 4 diagonal un- 



After these equations are solved for the diagonal elements, 
the unknowns of group (c) are obtained easily. 

This method becomes an effective tool for the solution 
of the matrix equation for systems of any order when its fol- 
lowing general properties are realized: 

T. The simplicity of the operations in the derivation of 
the final set of n equations is the same for systems 
of any order. These operations are simple substitu- 
tions of one element by another, which is done by 
manipulation of the unknowns indices. Relatively few 
arithmetic operations are done in arranging the final 
set of equations. 

2 . Because of the special arrangement of the unity ele- 
ments in. the matrices A and A T , there are very simple 
relations between the indices of the unknowns in each 
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equation, which enable one to describe the reduction, 
operation that was discussed above for the 4 ' th order 
case in the following four steps general scheme: 


The encircled elements are 

X 11 © *13 x 14 x 15 X 16 

the elements of group (a) which 

x 22 ^23) X 24 X 25 X 26 

s — s. 

are equal to corresponding ele- 

X 33 X 35 X 36 

ments of the matrix X n divided 

x 44 (*4| X 46 

by 2. 

X 5 5 © 



X 66 


b) Consequently, the group (b) 
elements are obtained. 



The group (c) elements are 

X 11 X 12,' X 13 X 14^ X 15 x 16 
* / 

expressed in terms of the 

x 22 X 23.- X 24 x 2 5,^26 

diagonal elements. 

x 33 X 34,^35 X 36 


X 44 X 45^ X 46 


X 55 x 56 

• 

x 66 


d) All elements are substituted into the group (d) equations 
that become the reduced set of n equations for the solution 
of the n diagonal elements. 
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In the machine procedure steps 3 and 4 will be built by index 
manipulation in the instructions that form the reduced equa- 
tions set. 

After the reduction it is left to solve n linear equation 
in n unknowns, or, in other words, it is left to solve an 
equation of the type: 

B x = c {4.1} 

where x and c are n dimensional vectors and 3 is an n x n di- 
mensional matrix, This reduction can be regarded as a reduc- 
tion of the matrix equation (1.6) to the matrix equation (4.1), 
or as a reduction of a set of n x n. linear equations to a set 
of n linear equations. For comparison, in the. method suggest- 
ed by Palsson the equation: 

E x, = } E. ,x n (2.3) 

n-1 , 1 - 1—0 

x=l n-i 

t ' r •/ *’ 

(which can be regarded as another reduced form of the ori- 
ginal matrix equation) must be solved for the first column of 
the state matrix X. Numerous arithmetic computations must be 
done to obtain the E^ matrices first, and then the sum vector 
on the right hand side of the equation (see sections 2.2 and 
2.3). The number of operations that must be done to obtain 
the matrix B and the vector c in equation (4.1) is rela- 
tively small. There is no loss of accuracy at all in obtain- 
ing the matrix B, and there is a minimal amount of arithmetic 
computation in obtaining the vector c. 
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If the operations that must be done prior to the solu- 
tion of equation (2.3) are compared to the reduction opera- 
tions that must be done prior to the solution of equation 
(4.1) in the new method, the advantage of this new method 
can be realized. The computation of the remaining elements 
of the X matrix which is done by simple substitutions in the 
new method, is also much more complicated in Palsson’s method, 
which again involves matrix operations. (as was shown in 
section 2.5, if the Model Performance Index is used and the 
order of the model is small enough, the computation of the 
remaining elements of the X matrix can be done by simple 
shifting operations, in which case this part of the computa- 
tion in the two methods will be of comparable simplicity.) 

4 . 3 An Alternative Approach to the Derivation of a Linear 

System Transfer Function 

One method for the derivation of the transfer function 
of a multi-loop system with a multi-input multi-output con- 
trolled member was mentioned in section 1.4. For high order 
systems the matrix operations of equation (1.21) require a 
great deal of computation adjoined by loss of accuracy. As 
was discussed before, significantly inaccurate system repre- 
sentation may not be tolerated by the optimization procedure. 

The approach suggested in this section is to obtain the 

transfer function of the controlled member first (if it is 

not given in this form) and then obtain the transfer function 
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of the whole system by multiplication of the system members 
transfer functions in the proper order. This approach 
actually replaces the matrix operations by polynomial op- 
erations. The accuracy of arithmetic operations between poly- 
nomials has not been examined in this work, and while the 
amount of computation would probably be smaller than in the 
n dimensional matrix operations, improved accuracy is not guar- 
anteed. When the system members are represented by their 
transfer functions, numerator and denominator polynomials 
can be multiplied and summed in any correct order, to give 
the over-all transfer function (e.g. obtaining the transfer 
function for each loop, proceeding from inner to outer loops). 

If the controlled member is given by a set of dynamic 
equations, its transfer function must be derived first. In 
the case of a single input this can be done by application of 
equations (1.20) and (1.21). (Notice that in many practical 
cases the order of the controlled member may be appreciably 
smaller than the order of the whole system, and the small di- 
mensional matrix operations to obtain its transfer function 
would not cause inaccuracies.) In the case of a multi-input 
controlled member whose inputs may be coupled, the transfer 
function can be derived by polynomial manipulation. Consider 
a controlled member with two inputs and three outputs, whose 
Laplace transformed dynamic equations are: 

AxYl + ^2 y 2 + -3 y 3 = Sl U]_ + JD 2 u 2 (4.2) 
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where y. is the i 1 th output and u. is the j 1 th input and 
and Ej are 3 dimensional vectors whose elements are poly- 
nomials of the Laplace variable. 

In the control system the controlled member outputs may 

be fed~bacl<~to“its — input s~th rough— con t-ro-l— components-. — All _ - 

the feedback paths between controlled member inputs and out- 
puts can be arranged so that there is not more than one trans- 
fer function in each path, as shown in the figure. 



Figure 4.1 A Controlled Member with Coupled Inputs 

If the transfer function between the input and any output 
is desired, then u^ is expressed as: 
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and substituted into equation (4.2): 


- 


G 12 )y l + 


^—2 G 22 5y 2 


(£3 


23 


)y 3 = ^ 1 U 1 


or 


— l y l + — 2 y 2 + — 3 y 3 ^1 U 1 


(4.4) 


k k k 

where the elements of the vectors £ 3 , and £3 are ratios of 
polynomials of the Laplace variable, since: 

^ = h - G 12 
-2 = -2 G 22 


£.3 - £3 - g 23 


(4.5 


The transfer function between the input u^ and any of the 
outputs, can be obtained by application of Cramer's rule, e.g. 


5l & 4*3 


u 


* * * 
-1 -2 -3 


(4.6) 


To obtain the transfer function as a ratio of two expanded 

polynomials, the numerator and denominator determinants must 

be expanded, which requires operations between the polynomial 

elements of the vectors m^ , £ ^ , £3 etc. In many practical 

cases some of the transfer functions G. . would be zeros or 

■L 3 

simple gains and the others would usually be ratios of low 
order polynomials, so that these operations would usually be 



simple and accurate. 

It is interesting to realize that the operations between 

the elements of the vectors actually replace the operations 

that would have to be done between the transfer functions 

G. . and the controlled member transfer function in the re- 
J 

gular case without coupled inputs, to obtain the over-all 
system transfer function. Thus, the coupling between the 
inputs does not require additional operations but rather a 
different arrangement of the controlled member equations and 
the corresponding control system members, according to the 
procedure that was described above. 

4 . 4 A Method for Finding the Relationships between the 

Transfer Function Coefficients and the Design Parameters 
In a parameter optimization procedure a representation 
of the' .control system must be done at each step of the opti- 
mum search in order to evaluate its performance and to com- 
pute the new step. The derivation of the system transfer 
function may be a time consumming process, especially for 
high order systems. In the algorithm that was discussed pre- 
viously the transfer function must be computed twice for each 
design parameter for the computation of the gradients at each 
step. Inaccurate computation of the transfer function co- 
efficients may cause problems like the ’’performance index in - 
crease in half a step” and incorrect gradient values, arising 

usually when the system transfer function is computed twice 

\ 
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for close values of design parameters. 

If the transfer function coefficients can be expressed 

\ 

as simple functions of the design parameters , then after 
these functions have been determined, the values of the co- 
efficients at each step of the optimum search can be found 
with relatively small additional ‘effort, and the- derivatives _ 
of the system coefficients with respect to the design para- 
meters can be found by directly deriving these functions. 

The method proposed in this section applies to high 
order complicated systems , where the computation of the 
transfer function must be done numerically and the distribu- 
tion of the design parameters cannot be followed easily, if 
at all, through the computation from the initial system, 
given by its components, to the final mathematical representa- 
tion. 


The design parameters may be system gains , time constants , 
damping ratios and natural frequencies. The over-all system 


transfer function will be considered in the following form: 


G(S) = 


, ~m . , 0 m-l . , 
b S + b •, S + . . . +b n 
m m-1 0 

a S + a *^+ . . . +a n 

n n-1 0 


The coefficients : b^ 0<i<m 

a . 0 < i <n 

can be replaced by the coefficients: c^ 0<k<m+n 

where: . ik=i 

c k = b i { 0 <i<m 

c = a (k=m+j+l 

k a j S0<j<n 
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When the system transfer function is derived by properly 
ordered multiplication of its members transfer functions, 
the over-all transfer function coefficients are sums of 
products of design parameters and constant numbers. Each 
one of the transfer function coefficients can be expressed 
by the following general expression: 


c. = s k.np. 

K V j 1 

3 J 


(4.7) 


where k. are constants and P. are design parameters. Ex— 

3 

pression (4.7) should be read as follows: "Each polynomial 
coefficient is equal to the sum of all the combinatorial pos- 
sibilities of products of the design parameters and constant 
numbers". If there are two design parameters, for example, 
each polynomial coefficient can be expressed as: 

c k = k l + k 2 P l + k 3 P 2 + k 4 P l P 2 (4 * 8) 

In the ; case of three design parameters the general expression 
of the coefficients would be: 


C k k l +k 2 P l +k 3 P 2 +k 4 P 3 +k 5 P l P 2 +k 6 P l P 3 +k 7 P 2 P 3 +k 8 P l P 2 P 3 ^ 4 * 9 ^ 


In the general case of t design parameters the number of 

products or the number of constants k^ can be computed by the 

combinatorial expressions : 

0 12 / 

C l + C l + ’ + c '= 2 

where 

r k - l? 

I ~ U-k) ! k! 
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Example 4.1. The transfer function of the system in 
Figure 4.3 is: 

G (S ) '= IQS + 20 

P 1 S 3 + S 2 + (P + 10P 2 )S + (1+20P 2 ) 

The constants k^ for each polynomial coefficient are listed 
in the following table: 


coef ficient 

k i 

CM 

k 3 

k 4 

c 0 =b 0 =20 

20 

0 

0 

0 

c l =b l sl ° 

10 

0 

0 

0 

C 2 =a 0 =1+20P 2 

o 

0 

20 

1 

C 3 =a l~ P l + 10 P 2 

0 

1 

O 

rH 

0 

C 4“ a 2“ 1 

1 

0 

0 

0 

C 5 =a 3 =P l 

0 

1 

0 

0 


Figure 4,2 The Constant k^ in the Polynomial Coefficients of 
_ Example 4.1 



Figure 4.3 A Control System for Example 4.1 
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To find the relationship between the design parameters 
and the transfer function coefficients it is necessary to 
find the values of the constants of the general equation 

(4.7) fo r each polynomial coefficient. Notice that the equation 

(4.7) is linear in the constants (also see equations (4.3) 

and (4.9)). Thus, the solution for the constants of each de- 
sign parameter can be obtained by forming sets of enough 
equations for each transfer function coefficient, to solve 
for its constants. 

To solve the 2 n constants in a polynomial coefficient 
2 n equations are needed. These equations can be produced by 
substituting 2 n sets of design parameter values into the 
given control system, and computing the corresponding sets 
of polynomial coefficients. All the combinatorial products 
of each set of design parameters are also needed to form the 
equations. 

In the case of two design parameters, for example, four 
equations are needed to solve for the four constants in 
each polynomial coefficient. Four sets of design parameter 
values will produce the following set of equations for the 


constants in the k’th polynomial coefficient : 


= k + P k + p k- + P. 

(1) 1 X (l) 2 2 (1) 3 X (i) 2 (i) 4 


P„ k 
P„ k 
P. k 


k ^1 P 1 ^2 P 9 k 9 P -I *■ 9 a 

( 2 ) 1 ( 2 ) 2 2 ( 2 ) 3 X ( 2 ) 2 (2 ) 4 

= k l + P k 2 + P k + P * „ 

(3) 1 1 ( 3) 2 2 (3) 3 1 ( 3) 2 ( 3 ) 4 

= k + P k 2 + p 9 k ? + ?! P? k 4 

*(4) i (4) 2 2 ( 4 ) 3 1 ( 4 ) 2 ( 4 ) 4 


(4.10) 
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The solution of the constants can be obtained by Cramer's 
rule . e . g . : 


c k 

k (l) 

X (l) 

P 2 

Z (D 

P 1 P 2 
i (l) (X) 

c. 

p, 


P, p„ 

k (2» 

1 

(2) 

2 ( 2 ) 

X (2) 2 (2) 

c. 

P, 

P^ 

Pi P- 

k (3> 

1 ( 3 ) 

2 (3) 

1 2 
(3) ^ (3) 

C k 

K (4) 

P 1 

(4) 

P 2 

(4) 

P P 

1 ( 4 ) Z (4) 

k l~ 

D 

( 4 • i. X ; 

where : 




1 

P. 

P. 

Pi Pi 


1 (1) 

2 (1) 

hi) 2 (i) 

1 

P, 

P 0 

Pi Pi 


1 ( 2 ) 

2 (2) 

1 ( 2 ) 2 (2) 

D= 1 

P, 

P^ 

P P 


(3) 

(3) 

X (4) ^ (3) 

1 

P, 

P^ 

Pi Pi 


1 (4 ) 

2 (4 ) 

1 2 
1 (4 ) Z (4) 




(4.12) 

Notice that the divider D does not 

depend on the system at 

all, bgt only on the values 

of the 

design parameters that were 

chosen to produce the equations. Also notice that in the num- 

erator determinant 

only one 

column 

depends on the system, and 


the other elements are functions of the. chosen design para- 
meters. It follows that four sets of values for two design 
parameters can be chosen once and then be used for all systems 
with two design parameters. The divider can be computed for 
these constant values independently of the control system. 
Since the numerator determinants in the expressions of the 
constants differ only by the values in the column of system 

coefficients c v and by the location of this column, the co- 
K (i) 
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factors of all the elements of the determinant D can be 
computed as system independent constants and be used for 
the computation of each k^ . This computation will then be 
done by multiplying a set of system coefficients by the 
corresponding constant cofactors and dividing by the constant 
divider. This, of course, applies to any number of de- 
sign parameters as well. 

Example 4.2* Consider all systems with two design para- 
meters. Four sets, of design parameter values are necessary 
to produce the equations for the solution of the constants 
k j , as was explained previously. The following sets may be 
selected : 


set 

1: 

p i=i 

P 2 =1 

set 

2: 

p r 2 

< 1 

II 

Oi 

0* 

set 

3: 

i— J 

II 

V 2 

set 

4: 

P 1 =2 

p 2 =2 


The corresponding divider D (see equation 4.12) is 1 and the 
cofactors are listed in the following table: 


column 

cofactors 

1 

4 

2 

-2 

■ 

2 

2 

2 

-1 

B 

3 

-2 


2 

B 

4 

m 



B 


Figure 4.4 Cofactors for Systems with Two Design Parameters 
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Suppose the constants in the polynomial coefficients of 
the system in example (4.1) are to be determined by this 
method. The computation of the constants in the coefficients 
Cq and is illustrated below: 

c 0 = 20 

For the four sets of design parameter values we get: 






20 


Multiplying the coefficients column by the cofactors we get: 


k = 20M -20 *2 +20 (-2) '-20 (-1) = 20 

-k 2 = 20-2 -20*2 +20 (-1) -20 (-1) = 0 

k 3 = 2 0 (-2) -20 ( -1 ) +20 * 2 -20 *1 = 0 

-k 4 = 20 (-1) -20 ( -1 ) +2 0 * 1 -20 *1 = 0 


and for-: 


c 3 = p i + 10P 2 


the coefficients column as obtained by substituing the design 
parameters values into the control system and computing the 
transfer function coefficients is: 


c- = 1 + 10 = 11 
J (l) 

c, = 2 + 10 « 12 

J (2) 

c -j = 1 + 20 = 21 

(3) 

c = 2 + 20 = 22 

J <4) 
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0 


and the constants for are: 

k 1 = 11*4 - 12-2 + 12 (-2) - 22 (-1) = 

-k 2 = 11*2 - 12-2 + 21 { -1 ) - 22 (-1) =-l 

k 3 = 11 (- 2) -12 (- 1 ) + 21 • 1 - 22-1 = 10 

-k 4 - 1K-D-12 (~1)+21-1 - 22-1 = 0 

A program for the computation of c factors and dividers 
using specified sets of design parameters values is given in 
Appendix C. Values of cofactors and dividers are also listed. 
The program is written for systems of up to six design para- 
meters and the listed results are for up to four parameters. 
The computation was found to be time consuming as the number 
of design parameters gets larger. It should be noted, however , 
that the values of the cofactors and the divider must be ob- 
tained but once, and then be used for all systems of a given 
number /of design parameters. The program also punches out 
the values of cofactors, as listed in the table in Appendix 
C. For the specific choice of design parameters sets the 
divider value is 1 for any number of design parameters. To 
use these results for finding the constants in the transfer 
function coefficients according to the procedure described 
above, the same design parameters sets must be used to obtain 
the corresponding transfer function coefficients. The proce- 
dure listed in the main program may be used to obtain the same 
sets of design parameters. 

Since the values of cofactors and dividers have been ob- 
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tained independently of the control system itself {only the 
number of design parameters needs to be specified)'/ the 
system's mathematical representation in terms of the design 
parameters can be done by finding the corresponding numeri- 
cal values of the transfer function coefficients/ for the com- 
putations of the constants k^ in each coefficient. This 
means that the system transfer function must be derived num- 
erically a number of times, even when this proposed method 
is used. This, however, is done before the optimization 

process begins. For t design parameters the transfer function 

l 

derivation must be done 2 times. This number is small in the 
case of few design parameters (for three design parameters 
the transfer function must be computed eight times) . In the 
case of. many design parameters this number would get consid- 
erably large (32 for five design parameters) . It should be 
•f- 

note d , -how ever, that in the currently used method, in the 
case of five design parameters the transfer function must be 
computed at least 30 times in three optimization steps. The 
number of steps grows rapidly with the number of design para- 
meters. Notice that the number of steps also grows with the 
order of the system, and in the currently used method this 
means more computations of the system transfer function, 
while in the newly proposed method the number of times that 
the transfer function must be computed depends only on the 
number of design parameters. 


65 



An even more significant advantage of the new method 
would be the accuracy of the computation as outlined by 


the following three points: 

1* As can be seen in Appendix C for the specific choice 
of design parameter values (only values of 1 and 2 
are used) the divider is always 1 and all the cofactors 
are small integers so that the computation of the 
transfer function coefficients involves minimal 
round-off errors imposed by the design parameters. 

This computation can always be checked for. accuracy 
by use of another set of design parameters values. 

The resulting kj constants must be the same for any 
choice of design parameters. If they vary for 
different sets of parameters values,' average values 
can be taken. 

2. Once the constants k^ have been determined for each 
tranfer function coefficient, the computation of the 
coefficients at each step involves very few arith- 
metic operations. Problems like the "performance index 
increase in half a step" are not likely to happen. 

■ 3. The computation of the derivatives of the transfer 

function coefficients with respect to the design para- 
meters is very simple. For example, in the case of 
three design parameters, using equation (4.9) we 


have : 
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k-> + k c P~ + 

2 5 2 6 3 

k-, + k-pT + k_P 0 

3 5 1 7 3 

k 4 + k 6 P l + k 7 P 2 


+ k g P 2 P 3 


+ k 8 p 1 p 3 


+ k 8 P l P 2 


Finally, it should be noted that if the system transfer function 
is computed by the state-space method (see equation 1.21), 
the transfer function coefficients discussed in this section 
appear divided by the coefficient a^. In this case the de-^ 
pendence of on the design parameters must be determined 
separately. If the polynomial approach is used for obtaining 
the transfer function (as proposed in the previous section) 
no special consideration must be given to a . 


67 


CHAPTER 5 


CONCLUSIONS AND RECOMMENDATIONS 


5 . 1 Conclusions 

The numerical accuracy problems arising in some stages 
of a parameter optimization technique have been presented 
and analysed. Special attention was given to the solution 
of the matrix equation ( 1 . 6 } and to the system mathematical 
representation. The accuracy obtained at each stage of the 

A' 

computation depends on the operations and the numerical 
values involved. 

A slight difference between the computed static sensi- 
tivities of the control system and the model, when the Model 
Performance Index is used, was found to cause some of the 
severe inaccuracy problems encountered in the design by 
parameter optimization. In the example case ’the difference 

7 

was one part of 10 . Another major source of inaccuracy 

is the inversion of the matrix E . Large values of the matrix 

n 

elements, which result in a large determinant value, may 
cause severe inaccuracies even when digital overflow or under- 
flow do not occur in the matrix inversion. This problem can 
be handled by reduction of the matrix elements before its 
inversion. The required reduction factor would depend on the 
transfer function coefficients and on the order of the system. 
No specific criterion for the value of the reduction factor 
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has been established. However, keeping the matrix determin- 
ant value below 10^ has given satisfactory results in the 
cases that were examined. A situation where the diagonal 
elements are significantly smaller than elements on their 
right in the rows is also unfavorable. Scaling the system 
coefficients is a useful means for conditioning matrices 
and vectors in a desirable way. The importance of accurate 
system representation to the optimization procedure has 
been emphasized. 

Alternative approaches to the problems of the system 
mathematical representation and the solution of the matrix 
equation (1.6) have been suggested for better accuracy and 
cost effectiveness. 

5.2 Recommendations 

Some of the changes recommended here have been applied 
to the parameter optimization program used at the Measure- 
ment Systems Laboratory at M.I.T. and have given the desired 
improvement in the accuracy of the computation. The recom- 
mendations are listed below: 

1. Fixing the system static sensitivity to be exactly 
equal to that of the model (fixing the value of 

°i 

at -1 when the model static sensitivity is 1) has 
resulted in immediate improvement of the computed 
results . 

2. Reducing the elements of the matrix E^ prior to its 
inversion when these elements are large has improved 
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the accuracy appreciably. The decision to reduce the 
matrix elements and the choice of the reduction fac- 
tor can be built in the program by examination of the 
characteristic coefficients and the system order. 
Consideration of the general expression of the E n 
matrix (see Appendix A) may be useful in estimating 
the magnitude order of the matrix determinant 
without computing the determinant value itself. 

3. Scaling the system coefficients can be used to con- 
trol the conditioning of matrices and vectors. It 
is recommended that the scaling operation be sepa- 
rated from the reduction of the matrix elements that 
was discussed above. Additional research should 
establish more specific criteria for scaling the 
system coefficients on a selective basis. Condi- 
tioning the matrix E^ prior to its inversion is one 
use of such scaling. The reduction of the matrix 
elements, on the other hand, should be done routinely 
whenever these elements are too large. 

4. The option of inverting the matrix in extended 
precision has been suggested. The matrix inversion 
program and the control cards are given in Appendix B. 

5. The computation of the initial conditions can be done 
in double precision without substantial additional 
storage. This change has given significant improve- 
ment in the computed results. 
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6 . As was dis cussed in section 2 . 5 some of the computa- 
tion currently done in the determination of the ele- 
ments of the state matrix X is unnecessary, since 
many of these elements do not take part in the com- 
putation of the performance index . Except for the 
first column of the matrix X only the elements in- 
cluded in the (£+l)x(£+l) left upper corner of the 

X matrix must be computed (t being the order of the 
model) * 

7. In the optimization program currently used a routine 
for checking and correcting the matrix X to be sym- 
metric is include* This is not neccesary since by 
the shifting operation in the computation of the 
columns of the matrix, discussed in section 2.5, 

the computed X matrix is always symmetric. 

8. The method for solving the matrix equation (1.6) , 
that was proposed in section 4.1 would provide a more 
accurate and less expensive solution than the method 
currently used. This method can replace the current 
algorithm without additional changes in the optimiza- 
tion program. 

9. The method proposed in section 4.2 for the computation 
of the system transfer function may be more accurate 
than the one currently used, but its accuracy has 

not been examined. It suggests a different approach, 
avoiding high order matrix operations. 
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The method proposed in section 4.4 is recommended 
as an accurate and cost effective way for comput- 
ing the system transfer function and the gradients 
at each step of the optimization procedure. 



APPENDIX A 


THE MATRICES E i 


The first three matrices are given below: 
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W S3 



a ~ 
n- 3 

-a « 
n-2 

a n-l 

-1 

... o 

0 . . 

. .0 

0 

n- J 

" a n-2 

a n-l 

-1 

0 

0 

0 

0 

a n-3 

" a n-2 

a , 
n-1 

~ 1 

0 


0 

0 

a , 
n-1 

-1 

a o 

a l 

0 

2a n-l 

-2a ^ a 

a -2a. a 


-2a 2 

0 n-1 

c 

h- 

y 

1 

H 

n-3 n-1 n-2 

n-1 

2 a A a 1 
0 n-1 

-2a~a .. +2a n a 2 , 

0 n-1 1 n-1 

+2a n-2 a n-l 

2a n-3 

2 ~n-2 a n- J 
+ 2 

a i 
n-1 


The first and last columns of and are given below: 


a n-4 

0 

« 

0 

• 

1 

0 

1 

“ a o 

-2a , 

n-1 

2a 0 a n-l 

2a n-l 

2a 0 a n-l 

~ 2a n-3 +2a n-2 a n-l~ 2a n-l 

3 

2 3.a 3 “2 9, ~ -j 2 S u *i * • * 

0 n-3 0 n-2 n-1 0 n-1 

••• 4a n-l a n-3- 4a n-2 a n-l +2a n 
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As the computation of the matrices proceeds , the elements 
shift upwards in the columns, inverting signs, while the 
last element of, say, the }<!th column is computed by multi- 
plying the k'th column of the previous E matrix by the last 
row of the system's matrix A. 
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APPENDIX B 


A PROGRAM FOR 

MATRIX INVERSION IN QUADRAUPLE PRECISION 


The program is a modification of the IBM/I Scienti- 
fic Subroutine Package MINV program for extended precision 
and for utilization by a FORTRAN main program. Also 
listed is the control cards set up for the IBM/370 system 
at M.I.T. Information Processing Center. This set up may 
change due to improvements of the system. 
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3 . X The Matrix Inversion Program 


MINV* • MINV 10 
/*****#*••** ##*******#*##*#**##**** *<mhmmhhmhhhhhhhhhhhh»**<imhmhmm>*****/M I NV 20 

/* */MINV 30 
/* TO INVERT A MATRIX */MINV 40 
/* */MINV 50 




PROCEDURE (AO.N*DO*CONQ) OPTIONS (FORTRAN) *• 
PUT SKIP LIST (* JUST ENTERED MINVm 
PUT SKIP DATA ( AQ *N*DQ *CONQ ) 1 

*** 

7 

DECLARE 

MINV 

80 

ERROR EXTERNAL CHARACTER ( 1 ) * 

MINV 

90 

(I*J*K*N*L(N) *M(N) ) 

/* 

MINV 

100 


FIXED BINARY* 


MINV 110 


*/ 

FIXED BINARY (31*0) * 

(BIGA*HOLO*D*CON*S) BINARY FLOAT(109>* 

A(UN»1SNJ BINARY FLOAT (109) * 

OC ( 1 »N* 1 :N) BINARY FLOAT (109) * 

(DQ*AQ (20*20) *CONQ) BINARY FLOAT (53) *. /^DOUBLE PRECISION INPUT ***131 

/* BINARY FLOAT (S3)** /*DOUBLE PRECISION VERSION /*0*/MINV 140 

/* */MlNV 150 

D=DQ * • 

00 1=1 TO N* • 

DO J=1 TO N* * 

A < I * J) =AQ < I * J) *. 

END** 

' END** ' ."v ; v 

CON=CONQ*. 


ERROR=* 0 f * • 


MINV 160 

IF N LE 0 


MINV 170 

THEN DO** 


MINV 180 

ERROR= • 1 * * • 

/* ORDER OF MATRIX « 0. 

*/MINV 190 

GO TO FIN*. 


MINV 200 

end*. 


MINV 210 

IF CON= 0 


MINV 220 



-Jf 

CO 


THEN 
/*THEN 
ELSE 
IF N 
THEN 


= I • OE-5 , • 
-1 ♦ OE-1 5 » * 
=CON t . 


= I/O 


» • 


N 


» • 


s 
s 
s 

= 1 
DO,. 

0 =A ( 1 , 1 ) , . 

IF ABS(O) LE S 
THEN 00,. 

ERR0R=*2« , 

ENO, • 

ELSE A (1,1) 

GO TO FIN,. 

END,, 

= 1 . 0 *. 

00 K = 1 TO 
L (K ) =K , * 

H (K ) =K , » 

BIGA =A(K,K) 

DO I=K TO N*. 

00 J=K TO 
IF ABSOIGA) 

THEN DO,. 

8IGA =A(I,J) ,. 
L(K) =1,. 

M(K) =J,. 

END,, 

ENO*. f 

END*. 

J S L (K ) , . 

IF L <K) GT K 
THEN DO,. 

DO I * 1 TO N,. 

HOLD *-A(K,I)t. 
A(K,I)=A(J,I) , . 
A(J,I)=HOLO,. 

END,. 

END*. 

1 =M (K ) , , 


/* SEARCH FOR LARGEST ELEMENT 


N, 


LT A8S ( A ( I , J) ) 


/* single precision VERSION /*S*/MINV 236* 
/* DOUBLE PRECISION VERSION /*D*/MINV 240 

MINV 250 

/* INVERT A SCALAR */MINV 260 

MINV 270 
MINV 280 
MINV 290 
MINV 300 
MINV 310 
MINV 320 
MINV 330 
MINV 340 
MINV 350 
*/MINV 360 
HINV 370 
MINV 380 
MINV 390 
MINV 400 
MINV 410 
MINV 420 
MINV 430 
MINV 440 
MINV 450 
MINV 460 
MINV 470 
MINV 480 
MINV 490 
MINV 500 
*/MINV 510 
MINV 520 
MINV 530 
MINV 540 
MINV 550 
MINV 560 
MINV 570 
MINV 580 
MINV 590 
*/MINV 600 


/* INTERCHANGE ROWS 


/* INTERCHANGE COLUMNS 



/♦ 

/♦ 

/* 


IF M (K) GT K 
THEN DO*. 

DO J = 1 TO N». 

HOLO =-A ( J*K) v, 

A { J*K ) = A ( J * I ) *. 

A< J*I)=HOLD*. 

END* • 

END*. 

IF ABS(BIGA) LE S 
THEN DO*. 

D =0*0** 

GO TO COMP,* 

END*. 

DIVIDE COLUMNS BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS 
CONTAINED IN BIGA) 

DO I = I TO N*. 

IF I NE K 

THEN A(I,K)*A(l'»K)/(*A(K*Kj ) *• 

END** 

DO I a I TO N*. /* REDUCE MATRIX 

IF I NE K 
THEN DO*. 

DO J a 1 TO N*. 

IF J NE K 

THEN A(I*J)=AU*K)*A(K*J)*A(I*J) 

END*. 

END*. - 

END*. 

DO J s I TO N*. 

IF J NE K /♦ DIVIDE BY ROW PIVOT 

THEN A (K * J) cA (K*J)/A(K»K> * . 

END*. 

D =D*A<K*K>*. /* COMPUTE DETERMINANT 

IF ABS(D) LE S 


MINV 610 
MINV 620 
MINV 630 
MINV 640 
MINV 650 
MINV 660 
MINV 670 
MINV 680 
MINV 690 
MINV TOO 
MINV 710 
MINV 720 
MINV 730 
♦/MINV 740 
♦/MINV 750 
♦/MINV 760 
MINV 770 
MINV 780 
MINV 790 
MINV 800 
♦/MINV 810 
MINV 820 
MINV 830 
MINV 840 
MINV 850 
MINV 860 
MINV 870 
MINV 880 
MINV 890 
MINV 900 
♦/MINV 910 
MINV 920 
MINV 930 
♦/MINV 940 
MINV 950 
MINV 960 


COMP.. 



FINAL ROW AND COLUMN INTERCHANGE 


00 

o 


THEN DO*. 

ERROR- *2* *• 

GO TO FIN*. 

END*. 

A(K*K)=1.0/A(K,K) *. 

END*. 

/* 

✓* 

/*• 

K =N * • 

LOOP.. 

K =K - 1 * . 

IF K GT 0 
THEN DO*. 

I -LOO*. 

IF I GT K 
THEN DO*. 

DO J =1 TO N*. 
HOLD =A(J*K> *. 

A ( J*K) =-A ( J* I ) ♦. 
A < J* I) =HOLD * . 
END*. 

ENO*. 

J =M(K>*. 

IF J GT K 
THEN DO*. 

DO I = 1 TO N*i 
HOLD =A(K*I) ♦. 
A(K*I)=-A(J*D*. 
A(J*I)=HOLD* . 
END*. 

ENO*. 

GO TO LOOP*. 

END*. 

FIN.. 


/* DETERMINANT IS ZERO 


/* 


MINV 970 
*/MINV 980 
MINV 990 
MINV 1000 

REPLACE PIVOT BY RECIPROCAL */MINV10l0 

MINV1020 
»/MINV1030 
*/MINV1040 
*/MINV1050 
MINV 1 060 
MINV 1 070 
MINV1080 
MINV1090 
MINV 1100 
M I NV 1110 
MINV1120 
MINV 1130 
MINV 11 40 
MINV1150 
MINV 1160 
MINV 11 70 
MINV1180 
MINV 1190 
MINV1200 
MINV 121 0 
MINV1220 
MINV1230 
MINV1240 
MINV 1250 
MINV 1260 
MINV1270 
MINV1280 
MINV 1290 
MINV 1300 
MINV1310 


OQ =0 * • 



AQ=A * « 

■*/ 

DO I - I TO N* 

DO J = 1 TO NJ 
AQ ( I * J) = A ( 1 1 J > ; 

end; 

end; 

PUT SKIP LI5TPABOUT TO RETURN FROM MINVM? 

CLOSE FILE (SYSPRINT) t 

RETURN** ' MINV1320 

/*ENO OF PROCEDURE MINV */WINV1330 


CO 



B,2 Control Cards Set Up for IBM/370 System at M.I.T. 


// »NAME • *CLASS=A*REGION=350K 

/miTID USER- 

/*MAIN TIME-2 * LI NES = 6 *CAROS = 0 
/./SXEPl~~.E X E ~ 

//C.SYSIM OD * 


t 

FORTRAN cards 

i ■ 


—//STEP 2 — EXEC- PLIKCL6*PARM # C=»CS (48 )* OBJECT ♦ NORUN ? COMPATIBLE * • 
// PARM.6=»C0MPATIBLE* 

y/C , S YS IN — DO — *,-&C a- BLKS-I-Z E *2 OO-O— - -2- 


t 


J?L/I cards 


i 


//L.SYSLIB DO DSN=U*M8230 • 7162* SU8R» PALSS * D I $P=SHR 
yy— oq ____ 


// 

y/- 


-D D-D S N=SY32.S S P- . S U B R ♦ Q IS P = S H R_ - 
DO DSN = SYS1.F0RTLIB*0ISP = SHR , H 


//{_ • SYS IN DO * 


^SUBfUDISP-SHR 


t 

OBJECT cards 


ENTRY MAIN 

NAME USER PRQG 

//G.FT06F00I DD SYSOUT = A*DCB={RECFM = VBA*LRECL=I37« 8LKS I ZE =2036 ) 
//G»FT05FQQJi — DQ »_t.D£B = BL K S 1 2 E= £0.0.0 _ 


! 

DATA cards 
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APPENDIX C 


THE SYSTEM INVARIANT DETERMINANTS AND COFACTORS 

A program for the computation of the system invariant 
determinant and its cofactors for the determination of the 
relationship between the system coefficients and the design 
parameters is listed. The program computs determinants and 
cofactors for systems of up to six design parameters. Also 
listed are the results for specified sets of design para- 
meters, that can be used for the computation of the constants 
k j , as described in section 4.4. ... 
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C.l The Computer Program 


DIMENSION P ( 6 ) * A ( 64 ♦ 64 ) * B < 64 * 64 ) * DET ( 64 ) 

1002 FORMAT ( * * ♦ ’ DETERMINANT =**F8*1/* THE COFACTORS OF THE ! COLUMN 
10F THE NUMERATOR DETERMINANT**/' CORRESPONDING TO THE CONSTANT K(I 
2)* ARE LISTED IN THE I GROUP*//* I*/) 

1003 FORMAT ( • 1 * * • NUMBER OF DESIGN PARAMETERS = **I1/* THE FOLLOWING SET 
IS OF DESIGN PARAMETERS VALUES ARE USED* ) 

1004 FORMAT (I3*3X* 10(F8.1*2X) / (6X * 10 (Ffl. 1 *2X> ) ) 

1005 FORMAT (10F8.1) 

1006 FORMAT ( * ***SET **I2** =**6<F7.1>) 

00 410 L = 1 * 6 

WR I TE (6* 1 003) L 
Ll=L-l 
L2=L-2 
L3=L-3 
L4=L-4 
11 = 1 

DO 10 I 1=1 *L 
10 P(I1)-1. 

WRITE (6* 1006) II* (P(J) *J=1*L) 

CALL AMAT (P*II*A*L) 

DO 20 11 = 1 * L 
1 1 = 1 1 ♦ 1 
P (1 1 > =2# 

WRITE (6*1006) II * IP <J) * J=1*L) 

CALL AMAT (P* I I *A*L) 

20 CONTINUE 
" IF(L.EQ*1) GOTO 200 
DO 30 I1=1*L1 

I 11=11 *1 V 

DO 30 12=1 1 1 *L 

II=II*1 

P(I1>=2. 

P < 12) =2# 

WRIT£(6*1006) II* (P(J) * J=1 *L) 

CALL AMAT <P* 1 1 * A*L) 

30 CONTINUE 




oo 

tn 


IF <L*EG*2) GOTO 200 
DO 40 1 1 = 1 , L 2 
111 = 11+1 
DO 40 12=111, LI 
121 = 12 + 1 
DO 40 13=121, L 
11 = 11+1 
P ( 1 1 ) = 2 # 

P(I2)=2. 

P ( 13 ) =2 • 

WRITE (6, 1006) II, (P(J> , J=1,L) 
CALL AMAT < P , 1 1 , A ,L ) 

40 CONTINUE 

IF (L ,EQ, 3 ) GOTO 200 
DO 50 11=1, L3 
111 = 11+1 
DO 50 12=1 11, L2 
121 = 12+1 
DO 50 13=121, LI 
131=13+1 
DO 50,14=131 ,L 
11 = 11+1 
P(Ii)=2. 

P(I2)=2. 

P(I3>=2. 

P(I4)=2. 


WRITE (6,1006) II, tP(J) ,J=1,L) 
CALL AMAT (P, 1 1 , A,L) 

50 CONTINUE 

IF(L.EQ.4) GOTO 200 
DO 60 11=1, L4 
111 = 11+1 
DO 60 12=1 1 1 ,L3 
121 = 12+1 


DO 60 13=121, L2 
131 = 13+1 






DO 60 I4=I31*L1 
141=14*1 
DO 60 15=141 *L 
11 = 11*1 
P<I1)=2. 

P(I2)=2. 

P ( 13) =2. 

P(I4>=2. 

P < 15 ) =2* 

WRITE (6 ♦! 006) !!♦ (P(J) fJ=lfL) 
CALL AHAT (P» 1 1 ♦ AtL) 

60 CONTINUE 

IF(L.EQ.S) GOTO 200 
11 = 11*1 
'“P(U)=2* 

P(I2>=2. 

P ( 1 3> =2 • 

P ( 14) =2 • 

P < 15) =2* 

P { 16) =2* 

WRITE<6*1006> II* (P ( J) ,J=1*L> 
CALL AMAT (P* II *A*L) 

200 LL=2**L 
I = i 

DO 210 I 1=1 #LL 
DO 210 JJ=1 »LL 
210 0(11, JJ)=A(II*JJ) 

IX=0 

KC=LL 

J=0 

GOTO 300 

100 J=J*1 

IF(J.GT.LL) GOTO 410 
1 = 0 

101 1 = 1*1 
' IK=0 



DO 221 II*1*LL 
irilt.EQ.I) GOTO 221 
IK= IK ♦ 1 
JK=0 

DO 220 JJ-ULL 

IF ( JJ.EQ* J) GOTO 220 

JK=JK+i 

B(IK*JK)=A(II»JJ) 

220 CONTINUE 

221 CONTINUE 
KC=Lt-l 

300 IF (KC«EQ« 1 ) GO TO 31 . 

IREV=0 

DO 22 IT=1 »KC 
K~I T 

„ 9 IF<8(K*IT) >21 *11*21 

o n i 

IF (K-KC) 9*9*51 

21 IF { I T-_K ) 12*14*51 

12 DO 13 M-l.KC 
TEMP=8 ( IT *M) 

B<IT*M)=8(K*MJ 

13 B<K*M)=TEMP 
IREV=I9EV*1 

14 IS=IT*i 

IF ( IS.GT *KC) GOTO ZZ 
DO 17 M=IS *KC v 

18 IF t B ( M * I T ) ) 19*17*19 

19 TEMP=B (M*IT)/B{IT*IT) 

" DO 16 N-IT t KC 

16 B(M*N)=B(M*N)-BUT*N)*TEMP 

17 CONTINUE 

22 CONTINUE 
DET ( I ) =1 • 

DO 2 I T=1 *KC ... ■- 

- 2 * DET ( I ) =DET ( I > *B ( I T ♦ IT ) 


OET ( I > = (-1 *00) **IREV*DET ( I ) 
GOTO 310 
51 DET ( I )*.0« 

GOTO 310 
31 OET ( I ) =B ( 1 * 1 ) 

310 IF<IX,EQ.0> GOTO 400 
' IF < I • EQ*LL) GOTO 320 
GOTO 101 

320 WR I TE ( 6 * 1 0 04 ) J * ( DET ( I ) * I = 1 1 LL ) 
WRITE (7,1005) (OET(I) ,I=l#LL) 
GOTO 100 

400 WRITE(6f 1002) DET < 1 > 

IX = 1 

GOTO 100 
410 CONTINUE 
STOP 
END 

SUBROUTINE AMAT (P, II ,A*L) 
DIMENSION A (64 164) , P (6) 

L1=L-1 

L2=L-2 

L3-L-3 

L4=L-4 

JJ=1 

A(il*JJ)*l. 

DO 10 1 1 = 1 *L 
JJ=JJ*1 ... *? 

10 A ( 1 1 ♦ JJ) =P ( 1 1 ) 

IF(L.EQ.l) GOTO 100 
DO 20. 11=1* LI 
111=11*1 
DO 20 12=1 1 1 * L 
JJ=JJ*1 

A(II*JJ)=P(I1)«P(I2) 

-?0 CONTINUE 



IF(L.EQ.2> GOTO 100 
00 30 1 1 -1 ,L2 

m = n*i 
do 30 i 2 =iu,li 
121 - 12+1 
00 30 13=121, L 
JJ=JJ*1 

A(II,JJ>=PCI1)*P(I2>*P(I3> 

30 CONTINUE 

IF (L «EQ« 3 ) GOTO 100 
DO 40 11=1, L3 
111 = 11*1 
DO 40 12=111, L2 
121 = 12*1 
00 40 13=121, LI 
131=13*1 
DO 40 14=131, t 
JJ = JJ *1 

A<n»JJ>>P<Il>*P<I2>*P(I3>«P(l4) 

40 CONTINUE 

IF (L.EQ.4) GOTO 100 
DO 50 I1=1,L4 
111 = 11*1 
DO 50 12=111, L3 
121 = 12*1 
DO 50 13=121 ,L2 
131 = 13*1 f- 

DO 50 14=131, LI 

141 = 14*1 HU- 

00 50 15=141, L * 

JJ=JJ*l 

* A(lItJJ)=Plll)*P(I2)*P(l3)*P(I4)*P<I5) 
IFIL.EQ.5) GOTO 100 

A(II»JJ>=P(1>*P<2)*P<3)*P<4)*P<5>*P<$) 
50 CONTINUE 
100 DO 500 1 = l ,L 



500 P<I>=1. 

RETURN 
- • • END 


CD 

O 



C.2 


Determinant and Cofactor Values for Specified Design Parameters 


NUMBER CF CESIGN iPARAMFTERS = 1 

THF FCLICMNG SETS CF CESIGN PARAMETERS VALUES A PE USEC 
SET 1 = l.C 

SEI 2 » i.O 

DETERMINANT* 1-0 

THE COFACTOaS CF THE I CCLUMN CF THE NUMERATOR CETERM IN/NT ■ 
CChRESPCNDING TO THE, CONSTANT Mil, ARE LISTEC IN THE I GRCLP 

I 

1 2.C 1.0 

2 l.Q 1.0 


NUMEER CF 

CESIGN 

PARAMETERS = 2 



THE FCLLCMNG SETS U.F CESIGN P FRAME 

TEftS VALUES ARE USED 

SEI 1 => 

l.C 

1 .0 



SET 2 * 

2-C 

1.0 



SEI 3 * 

1 -0 

2.0 



SFT A 3 

2.0 

2.0 



DF TERM [ NAN T = 

1 .0 



THE CCFACTCRS OF 

THF, I CCLUMN OF 

THE NUMERATOR 

DETERMINANT . 

CCRPESPCNCING TC 

THE CONSTANT Kill* 

ARE LISTEC IN 

THE I CRCIP 

i 

1 

<r.O 

2.0 -2.0 

-l.C 


2 

2.C 

2.0 -1.0 

-1.0 


3 

-2.C 

-l.C ?.C 

l.'C 


4 

-1 .C 

-1.0 l.C . 

. l.C 

_ 



NUMBER OF CESlGN PAP AME TEKS = 3 

TEE FCLLCwlNG SETS OF DESIGN PARAMETERS VALLES ARE LSEO 


SET 

1 = 

1-0 1.0 

l .0 




SET 

2 = 

2.C 1.0 

1.0 

- . . - -- -■ 

- - - 

. — . 

SET 

3 * 

l.C 2.0 

l.C 




SET 

4 * 

1.0 1.0 

2.C 




SET 

5 = 

2. C 2.0 

l.C 




SET 

6 = 

2 . 0 ; l.o 

2 • C 




SET 

7 * 

1.0 2.0 

2.0 




SET 

P. = 

2.0 2.0 

2 .0 




DETERMINANT* 1-0 





Thfc 

CCFACTCRS OF TEE 1 

COLUMN OF 

THE NUMERATOR 

HE 1ERMIN AN T . 

CORRESPONDING TO TEE CONSTANT Ktl), 

ARE LISTED IN 

TEE I 

CRCLP 

i 

1 



E.C 4.0 

-4.0 

4. C 

2.C 

- 2 . C 

2 


4.C 4.0 

-2.0 

2.0 

2.0 

-2 • C 

3 


-4.C -2.0 

4. C 

-2. C 

-2.0 

l.C 

4 


4.0 2.0 

-2.0 

4 . C 

l.C 

-2 • C 

5 


2.0 ; 2.3 

-2.0 

1.0 

2.0 

-1.0 

6 


-2.C : -2-0 

l.C 

-2.C 

-l.C 

2.0. 

7 


2.0 1.0 

-2.C 

2 . C 

l.C 

— 1 ♦ c 

8 


l.C 1.0 

-1.0 

1.0 

1.0 

-l.C 







....... 




2.C 

1.0 

l.C 

1.0 

2.C 

-i.O 

2 . C 

1.0 

l.C 

l.C 

l.C 

-1.0 

2 . C 

1.0 

l.C 

1.0 






CO 

CO 


NUMPER CF DESIGN PARAMETERS = 4 

THE FULLCH1NG SETS OF DESIGN PARAMETERS VALUES ARE USED 


SET 

1 = 

l.C 

1.9 

l.C 

1.0 

SET 

2 * 

2.C 

1.0 

1 .c 

l.C 

St I 

3 * 

1.0 

2.0 

1 .0 

l.C 

SE1 

4 *= 

l.C 

1,0 

2 .0 

l .0 

SET 

5 * 

l.C 

1,0 

l.C 

2 . C 

SE I 

6 =t 

2.0 

2 .0 

1 .0 

l.C 

SE T 

7 * 

2.C 

1 .0 

2. C 

l.C 

SET 

8 = 

2.C 

1.0 

l.C 

2 . C 

S E T 

9 = 

1.0 

2.0 

2 .0 

l.C 

SET 

1C * 

1.0 

2.0 

L .0 

2.0 

SET 

11 * 

l.C 

1.0 

2. C 

P.C 

SET 

12 = 

2 .0 

2.0 

2.0 

l.C 

SET 

13 = 

2 .0 

2.0 

l .0 

2.C 

SET 

14 = 

2.C 

1 .0 

2.0 

2.0 

SET 

15 * 

1 .C 

2.0 

2 .C 

2 . C 

SE! 

16 = 

2.0 

2 .1 

2.0 

2.0 


DETERPI N A M = 1.0 x „ 

THE CCFACTCRS OF THE ' l CCLUVN OF THE NUMERATOR DETERMINANT » 
CCRRESPCND I NG TC THE CONSTANT K T I I t ARE LISTED IN THE I CRCl 


I 

1 

2 

3 

4 

5 

6 
T 

a 

9 

1C 

11 

12 

11 

14 

15 


16.0 

e .o 

-P.c 

8.C 

-P.C 

-4. C 

4.C 

2.0 

-2.0 

2 .0 

-2.0 

-l.C 

8 . C 

8.0 

-4. C 

4. C 

-4.C 

-4.0 

2.0 

2.0 

-2.C 

2 . C 

- l.C 

-l.C 

-8.C 

-4.0 

« .0 

-4 .0 

4.0 

4 ,C 

-2.C 

-2.0 

?. c 

.-l.C 

7.0 

1 .0 

8 .0 

4.0 

-4.C 

e. c 

-4.C 

-2.0 

4 . C 

2,0 

-1.0 

2 . C 

-2 .C 

-l.C 

-b . C 

-4.0 

4. 0 

-4,0 

8 . C 

2.0 

-4 .C 

-1.0 

2.C 

-2.C . 

2.C 

l.C 

-4.C 

— 4.0 

4.0 

-2.C 

P.C 

4. C 

-i.C 

-2.0 

2.0 

-1.0 

1.0 

l.C 

4 .C 

4.0 

-P.C 

4. C 

-2.C 

-2.C 

2.0 

2.0 

-i.c 

2.C 

-l.C 

-l.C 

-4 . C 

-4.0 

2.0 

-2.0 

4 .0 

P.C 

-2 .C 

-1.0 

2. C 

-P.C 

l.C 

l .0 

4 .0 

2.0 

— 4 . C 

.« 4 . C 

-2 .C 

-2.0 

2 . C 

2.0- 

-1.0 

ii.c 

-2.0 

-1.0 

-4 .C 

-2.0 

4. C 

-P.C 

4 .C 

2 .0 

-2 .0 

-1.0 

2 . C 

-l.C 

2.C 

l.C 

4.C 

2.0 

-2.0 

4.C 

-4.C 

-l.C 

4.C 

1.0 

-1.0 

2.0 

-2.0 

-1.0 

2 .C 

2.0 

-P.C 

2 ♦ C 

-l.C 

-2 .G 

1 .0 

2.0 

-1.0 

l.C 

-l.C 

-l.C 

-2.C 

-2.0 

2.0 

-1.0 

2.0 

2 « C 

-l.C 

- 1.0 

?. C 

-l.C 

l.C 

1.0 

2 .0 

2-0 

-l.C 

2.0 

-P.C 

-l.C 

2.C 

1.0 

-1.0 

2.0 

-1.0 

-1.0 

-2 -C 

-1.0 

2. C 

-2.C 

2.C 

1 .0 

-2.0 

-1.0 

1 . c 

-l.C 

2.C 

l.C 

-l.C 

-1.0 

1.0 

-l.C 

l.C 

1 * c 

-l.C 

-1.0 

[.Cl 

-l.C 

1 .0 

1 .0 


16 


4 . C 

-4.0 

4 . C 

-4.0 

-2.C 

2.0 

4 . C 

-2 .0 

-2.C 

4.0 

-z'.c 

2.C 

4 . C 

-2.0 

-2.C 

4.0 

2.C 

-l.C 

-L.C 

2.0 

2.C . 

-2.C 

2 . C 

-1.0 

-l.C 

P.C 

2.C 

. -2.0 

-l.C 

t.O 

-l.C 

l.C 


4.0 -4.0 

2.0 -2.C 

-4.C 4 . C 

4.0 -2.0 

-2.C 4 . C 

-2.0 2.C 

2.0 - 1 • C 

- 1 . C 2. C 

4.0 -2.0 

-2.C 4 . C 

2, C -2.C 

2.0 - 1.0 

-l.C 2 . C 

1.0 - 1.0 

-2.C 2.C 

-L.C l.C 
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